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An  optimized  variant  of  the  nested  Markov  chain  Monte  Carlo  [n(MC)2]  method  [J.  Chem.  Phys. 
130,  164104  (2009)]  is  applied  to  fluid  N2.  In  this  implementation  of  n(MC)2,  isothermal-isobaric 
( NPT )  ensemble  sampling  on  the  basis  of  a  pair  potential  (the  “reference”  system)  is  used  to 
enhance  the  efficiency  of  sampling  based  on  Perdew-Burke-Ernzerhof  density  functional  theory 
with  a  6-3 1G*  basis  set  (PBE/6-31G*,  the  “full”  system).  A  long  sequence  of  Monte  Carlo  steps 
taken  in  the  reference  system  is  converted  into  a  trial  step  taken  in  the  full  system;  for  a  good  choice 
of  reference  potential,  these  trial  steps  have  a  high  probability  of  acceptance.  Using  decorrelated 
samples  drawn  from  the  reference  distribution,  the  pressure  and  temperature  of  the  full  system  are 
varied  such  that  its  distribution  overlaps  maximally  with  that  of  the  reference  system.  Optimized 
pressures  and  temperatures  then  serve  as  input  parameters  for  n(MC)2  sampling  of  dense  fluid  N2 
over  a  wide  range  of  thermodynamic  conditions.  The  simulation  results  are  combined  to  construct 
the  Hugoniot  of  nitrogen  fluid,  yielding  predictions  in  excellent  agreement  with  experiment.  ©  2009 
American  Institute  of  Physics.  [DOI:  10.1063/1.3200904] 


I.  INTRODUCTION 

The  behavior  of  matter  subjected  to  extreme  conditions 
is  a  topic  of  abiding  interest,  having  direct  application  to 
shock  physics  as  well  as  the  planetary"  and  geosciences."  Its 
relevance  to  the  latter  has  assumed  new  importance  in  the 
wake  of  global  climate  change,  as  sequestration  of  C02  in 
underground  reservoirs  may  constitute  a  vital  component  in 
emission-free  processing  of  fossil  fuels.4  The  focus  here  will 
be  on  equilibrium  characterization  of  warm,  dense  fluids 
such  as  those  produced  by  detonation  of  high  explosive  (HE) 
compounds.  Specifically,  Monte  Carlo  (MC)  simulation"'6 
will  be  used  to  characterize  the  locus  of  accessible  shock 
states  (known  as  the  Hugoniot  locus7)  for  molecular  nitrogen 
fluid,  an  almost  ubiquitous  component  of  HE  detonation 
product  mixtures. 

Recently8  (hereafter,  we  will  refer  to  Ref.  8  as  Paper  1) 
we  reported  a  thermodynamic  optimization  procedure  de¬ 
signed  to  improve  the  sampling  efficiency  of  the  nested  Mar¬ 
kov  chain  Monte  Carlo  [n(MC)2]  method.9’10  The  n(MC)2 
procedure  partitions  a  standard  MC  simulation  into  a  refer¬ 
ence  system  defined  by  an  approximate  potential  and  a  full 
system  defined  by  an  alternative,  more  accurate  one.  A  se¬ 
quence  of  elementary  moves  (in  the  NPT  ensemble,  these 
correspond  to  single-particle  displacements  or  volume  ad¬ 
justments),  each  accepted  with  Boltzmann  weight  in  the  ref¬ 
erence  system,  is  used  to  build  a  many-particle  composite 
trial  step  taken  in  the  full  system.  A  slight  modification  of  the 
acceptance  criterion  for  these  extended  trial  steps  recovers 
Metropolis  sampling  of  the  full  potential  without  having  to 
evaluate  it  at  each  configuration.  The  approach  developed  in 
Paper  1  improved  the  efficiency  of  the  n(MC)2  technique  by 
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allowing  the  thermodynamic  states  of  the  reference  and  full 
systems  to  vary  independently  such  that  their  respective  dis¬ 
tributions  attained  maximal  overlap,  thereby  raising  the 
mean  acceptance  probability  for  trial  composite  steps.  Be¬ 
cause  the  purpose  of  Paper  1  was  merely  to  illustrate  the 
procedure,  the  reference  and  full  systems  were  defined  by 
model  potentials  not  differing  in  computational  expense.  The 
present  work  serves  both  to  extend  the  application  of  opti¬ 
mized  (o-)  n(MC)2  to  a  full  system  characterized  by  density 
functional  theory  (DFT)  (Ref.  11)  and  to  construct  the  shock 
Hugoniot  locus  of  fluid  N2  at  this  level  of  theory. 

The  following  section  describes  in  detail  the  potentials 
used,  an  analytical  pair  potential  for  the  reference  system  and 
DFT  for  the  full  system.  Section  III  reviews  the  principal 
expressions  making  up  the  o-n(MC)2  procedure,  referring 
the  reader  to  previous  works  for  their  full  justification.  Sec¬ 
tion  IV  outlines  the  continuum  theory  of  shock  waves  as 
expressed  in  the  Rankine-Hugoniot  relations  as  well  as  the 
implementation  of  these  relations  in  the  context  of  atomistic 
simulation;  this  linkage  is  then  illustrated  by  construction  of 
the  N2  Hugoniot  locus  from  the  o-n(MC)2  simulation  results 
discussed  in  Sec.  III.  The  final  section  summarizes  and  offers 
some  thoughts  on  future  implementation  of  the  procedure. 

II.  POTENTIAL  ENERGY  EVALUATION 
A.  Specification  of  the  full  potential 

DFT  replaces  the  antisymmetric  A'-electron  wavefunc- 
tion,  expressed  in  spin-spatial  coordinates  x, 

'P  =  ^(Xj.Xj,  ....X^v),  (1) 

with  the  spatial  density  of  a  single  electron, 
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p(r)  =  Nj  •■■J  |^(x1,x2,  ...  ,xN)\2ds1dx2--- dxN.  (2) 

Although  this  procedure  appears  to  entail  considerable  infor- 
mation  loss,  the  first  Hohenberg-Kohn  theorem  “  guarantees 
a  unique  map  from  the  density  to  the  ground  state  energy,  Eg, 

p  =>//=>  'P  =>  Eg .  (3) 

Eg  can  then  be  expressed  as  a  functional  of  the  density, 
Eg[p],  and  its  constituent  parts  decomposed13  in  the  follow¬ 
ing  manner: 

Eg  [p]  =  rs[p]  +  J\p ]  +  £Ne[p]  +  £xc[p]  ■  (4) 

The  various  contributions  to  the  total  electronic  energy  are, 
respectively,  the  kinetic  energy  for  a  gas  of  noninteracting 
electrons,  the  (mean-field)  Coulombic  repulsion,  the  Cou- 
lombic  electron-nuclear  attraction,  and  the  exchange- 
correlation  functional.  The  first  three  terms  are  analytically 
soluble  whereas  the  fourth  subsumes  all  quantum  effects  (in¬ 
cluding  real-time,  dynamic  electron  correlation,  the  Pauli 
principle,  the  self-interaction  correction,14  and  correlated 
electron  contributions  to  the  kinetic  energy)  and  has  so  far 
proven  strongly  resistant  to  universal  description.  Although 
exchange-correlation  functionals  constructed  entirely  from 
first  principles  have  been  reported,15  most  £xc[p]  are  built 
from  fits  to  a  set  of  experimental  and/or  ab  initio  data.  Dif¬ 
ferent  functionals,  then,  are  distinguished  primarily  by  their 
exchange-correlation  component,  £xc[p]-  The  second 
Hohenberg-Kohn  theorem  "  establishes  a  variational  prin¬ 
ciple  for  the  density  in  the  same  manner  as  for  the  wavefunc- 
tion,  and  the  Kohn-Sham  equations  ~  permit  iterative  solu¬ 
tion  for  a  set  of  one-electron  orbitals  in  the  same  manner  as 
the  Roothaan  equations16  in  Hartree-Fock. 

All  of  the  quantum  chemical  calculations  reported  below 
were  carried  out  using  GAUSSIAN03.  In  test  calculations  on 
a  periodic  cubic  cell  containing  100  N2  at  density  p 
=  1.46  g/cm3,  we  evaluated  the  performance  of  four  different 
Exc[p]:  BLYP,18’19  PBE,20  PW91,21  and  BVWN.18’22  All  pe¬ 
riodic  calculations  of  the  liquid-state  energy  were  performed 
at  the  r -point.2  ’  In  addition  to  the  relative  computational 
speed  at  the  point  of  basis  set  convergence,  we  compared  the 
computed  bond  length  and  frequency  of  the  isolated  mol¬ 
ecule  to  their  experimental  values. 

The  basis  sets  examined  cover  the  range  from  minimal 
[STO-3G  (Ref.  24)]  to  triple  zeta  with  polarization  ( d  and  /) 
functions,  the  latter  containing  70  basis  functions  for  each 
molecule.  There  is  no  evidence  of  ionization  at  the  densities 
of  interest  here,  so  diffuse  functions  were  neither  necessary 
nor  desirable;  energies  computed  using  the 
6-3llG(3df ,3pd)  basis"5  were  thus  treated  as  the  complete 
basis  set  limit  (CBSL)  for  each  functional.  Table  I  records 
the  percent  deviation  in  the  total  energy  (relative  to  that  of 
100  isolated  N2  molecules)  for  a  periodic  box  of  100  N2 
molecules  in  a  typical  liquid-state  configuration  drawn  from 
the  reference  distribution,  calculated  as  a  function  of  Exc[p] 
and  basis  set,  from  that  calculated  with  the  same  functional 
at  the  CBSL.  All  energies  were  corrected  for  basis  set  super¬ 
position  error  (BSSE),26  the  spurious  stabilization  conferred 
upon  an  atomic  center  by  basis  functions  “borrowed”  from 


TABLE  I.  Percent  deviation  of  the  BSSE-corrected  energy  of  100  N2  mol¬ 
ecules  at  density  of  p=1.46  g/cm3,  as  a  function  of  basis  set  and  /Zx(:[p], 
from  that  calculated  using  the  same  functional  at  the  CBSL  [defined  as 
6-311G(3d7,3prf)  (Ref.  25)].  Three  of  the  four  functionals  tested  are  con¬ 
verged  to  roughly  1%  using  the  6-31G*  basis.  Based  on  these  results,  the 
6-3 1G*  basis  set  was  used  in  all  subsequent  calculations  reported  in  this 
study. 


Basis  set 

BLYP“ 

PW9l“ 

PBEa 

BVWNa 

STO-3G“ 

-40.2 

-42.9 

-33.8 

-34.1 

6-3  lGb 

2.5 

4.6 

5.3 

0.8 

6-31G*a 

-1.1 

-3.1 

-1.4 

-1.4 

6-311G**c 

-1.0 

-2.0 

-1.1 

-0.7 

“See  text  for  reference. 
bReference  52. 
‘References  27  and  53 


another  atom.  Our  target  accuracy  for  the  nitrogen  Hugoniot 
was  —1%,  and  although  this  value  does  not  translate  directly 
to  an  accuracy  constraint  on  the  DFT  energy,  it  does  serve  as 
a  useful  guide.  All  four  functionals  tested  were  converged  to 
roughly  1%  at  the  6-3 1G*  (Refs.  27  and  28)  level,  although 
the  performance  of  PW91  is  noticeably  worse  than  that  of 
the  other  three.  Because  full  system  energies  computed  using 
any  of  the  four  functionals  paired  with  6-3 1G*  yielded  good 
accuracy  in  reasonable  time,  this  basis  set  was  singled  out  for 
further  tests  on  the  isolated  molecule. 

Ground  state  properties  recorded  in  Table  II  reflect  little 
variation  among  the  various  flavors  of  DFT,  although  it  is 
interesting  that  all  four  outperformed  MP2  (Ref.  29)  in  pre¬ 
dicting  l,  the  equilibrium  bond  length  (correcting  for  anhar- 
monicity  makes  the  MP2  value  even  worse  relative  to  those 
of  DFT).  The  rightmost  column  of  Table  II  compares  relative 
timings  for  each  functional  paired  with  the  6-3 1G*  basis  set, 
and  speedup  factors  were  defined  as  the  ratio  of  time  t 
needed  to  compute  the  energy  of  the  full  box  of  100  N2  to  a 
reference  value  rref, 


T 

speedup  =  — . 

Tref 


(5) 


The  reference  time  was  defined  to  be  that  required  by  the 
slowest  functional,  giving  it  a  speedup  factor  of  unity.  PW91 


TABLE  II.  Performance  of  several  exchange-correlation  functionals  in  com¬ 
parison  to  those  of  MP2  and  experiment.  Points  of  comparison  include 
ground  state  bond  length  /  and  (harmonic-approximation)  frequency  w  of  the 
isolated  N2  molecule,  as  well  as  the  relative  time  required  for  a  three- 
dimensional  periodic  single-point  calculation  on  the  fluid  sample  of  100  N2 
molecules.  All  calculations  used  the  6-3 1G*  basis  set.  Speedup  factors  are 
defined  in  the  text,  the  largest  value  corresponding  to  the  fastest  time  and  the 
slowest  assigned  a  factor  of  unity.  Based  on  the  results  of  this  simulation, 
the  PBE  functional  was  chosen  for  use  in  all  subsequent  calculations  re¬ 
ported. 


EXC 

1(A) 

oj  (cm  *) 

Speedup 

BLYP 

1.118 

2337 

1.00 

PBE 

1.117 

2360 

1.13 

PW91 

1.116 

2364 

1.14 

BVWN 

1.116 

2346 

1.09 

MP2 

1.131 

2175 

Expt“ 

1.094 

2359 

‘‘Reference  39. 
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was  fastest,  although  the  difference  with  PBE  was  practically 
negligible.  On  the  basis  of  its  being  well  converged  at  the 
6-3 1G*  level,  predicting  satisfactory  ground  state  properties, 
and  (almost)  recording  the  fastest  calculation  time,  the  PBE 
functional  was  chosen  for  all  subsequent  calculations. 

We  conclude  this  section  with  a  note  regarding  omission 
of  local  density  approximation  (LDA)  functionals.  Although 
not  shown,  we  also  tested  SLYP19'30  and  SVWN22'30  in  addi¬ 
tion  to  the  generalized  gradient  approximation  functionals 
examined  above.  Calculated  energies  dropped  steadily  on  ap¬ 
proach  to  the  CBSL,  actually  dipping  well  below  zero  in  the 
case  of  SLYP  A  negative  energy  would  imply  that  it  is  ther¬ 
modynamically  preferable  for  a  sample  of  nitrogen  gas  to 
compress  spontaneously  to  a  density  of  1.46  g/cm3,  clearly 
indicating  catastrophic  failure  of  LDA  in  describing  the 
dense  N2  fluid  sample. 

B.  Parametrization  of  the  reference  potential 

Vital  to  our  scheme  for  building  equilibrium  averages 
with  DFT  energies  was  a  rapidly  computable  potential  func¬ 
tion  capturing  enough  of  the  essential  physics  to  serve  as  an 


2  0-25 r2(/2  +  /2)  +  (-  1)>,7  •  1,.)  +  (-  1)V,7  •  1,.)  +  (-!)' 

rab  ~  2 

r0 


adequate  reference  for  the  DFT.  After  having  converged  an 
accurate  combination  of  exchange-correlation  functional  and 
basis  set,  we  parametrized  the  reference  potential  on  the  ba¬ 
sis  of  this  model  chemistry. 

The  exponential-6  (Buckingham)  potential  was  used  to 
describe  pairwise  interaction,  <pab,  of  atomic  sites  a  and  b. 
The  sites  reside  on  diatomic  molecules  i  and  j,  so  the  com¬ 
plete  interaction  for  a  pair  of  reference  system  molecules  is 
summarized  in 

<p{rab)  =  —(  6e“(1-^>-4)  (6) 

a~6\  rj 

and 

2  2 

<Pij  =22  <p{rab) .  (7) 

a=  1  b=  1 

The  site-site  radii  can  be  re-expressed  in  terms  of  the  center- 
of-mass  (c.m.)  separation  vector  and  bond  vectors  1,  and  1; 
(of  lengths  /,■  and  Z;), 


'0-5(1/  ■!;■) 


Bond  lengths  were  fixed  at  /=  1.117  A,  the  value  predicted 
by  PBE/6-31G*  for  the  isolated  N2  molecule.  Parameters  e, 
a,  and  r0  represent  the  depth  of  the  well,  the  steepness  of  the 
repulsive  wall,  and  the  position  of  the  minimum,  respec¬ 
tively;  their  fitted  values  and  the  process  used  to  obtain  them 
are  provided  below.  Although  Eq.  (6)  realistically  captures 
the  exponential  character  of  hard-core  repulsion  and  the  r~b 

L 


T 


FIG.  1 .  Diatomic  pair  configurations  used  to  parametrize  the  reference  po¬ 
tential  whose  behavior  is  illustrated  in  Fig.  2.  Bond  lengths  were  held  fixed 
at  /=  1.1 17  A,  while  c.m.  separations  rtj  were  swept  over  a  range  of  roughly 
1.6-6. 6  A,  depending  on  the  configuration  type  (see  text).  Configuration 
types  are  characterized  by  four  angles,  , \\ ,  02 where  6  and  x  repre¬ 
sent  rotation  of  molecules  1  and  2  within  and  out  of  the  plane  of  the  page. 
Clockwise  rotations  are  positive,  and  all  angle  values  are  zeroed  to  those  of 
the  P  configuration. 


a=l  V 


o=2 


b=  1 


b= 2 


0,0,*,*] 
2  2) 


P 


(0, 0,0,0) 


3 1 

dependence  of  dispersion  at  long  range,  it  also  diverges  to 
-oo  at  very  short  range.  Using  the  parameter  set  described 
below,  the  potential  function  defined  by  Eqs.  (6)-(8)  reaches 
a  maximum  value  of  9  X  104  K  at  rah~  1.15  A,  so  enforcing 
a  minimum  allowable  separation  rmin=  1.20  A  safely  avoided 
the  (unphysical)  strongly  attractive  forces  present  at  small 
values  of  rab. 

In  order  to  ensure  suitable  parameter  values,  Eq.  (6)  was 
fitted  to  a  set  of  PBE/ 6-3 1 G*  energies  for  the  N2  dimer 
relative  to  those  for  a  pair  of  isolated  N2  molecules.  The 
distribution  of  pair  configurations  sampled  in  a  fluid  was 
approximated  using  the  four  fiducial  configuration  types  de¬ 
picted  in  Fig.  1,  denoted  as  P,  L,  T,  and  X.  Center-of-mass 
distances  were  scanned  over  ranges  of  1.59-4.23  A  (A-), 
1.85-4.50  A  (P),  1.59-5.29  A  (T)  and  2.38-6.62  A  (L)  in 
roughly  0.25  A  intervals  and  the  difference  between  Eq.  (6) 
and  the  PBE  result  was  recorded.  The  weighted  sum  of  these 
differences,  squared,  was  then  minimized  over  the  domain  of 
the  parameter  set. 


4  N, 

min  Err=  mini  VU,/(‘P°FT(M)  -  <p“p~6(k,/))2 

a,s,r0  a,E,r0  ^  i=\  k=\ 

(9) 

where  /  runs  over  the  set  of  pair  configuration  types 
{P,L,T,X}  and  k  indexes  the  N )  discrete  values  of  riy  The 
weights, 
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FIG.  2.  The  pair  potential  defined  by  Eqs.  (6)-(8)  and  (12)  as  a  function  of 
c.m.  separation  for  the  diatomic  pair  configurations  illustrated  in  Fig.  1 
(X,  T,  L,  and  P)  and  a  fixed  bond  length  /=  1.1 17  A.  Comparison  of  these 
curves  with  the  purely  spherical  potential  (/=0  A)  reveals  the  anisotropic 
character  of  (ptj.  Note  that  the  ordinate  is  plotted  on  a  log  scale. 


1 

WU  (max{<p?FT(k,/),9sw})2’ 


(10) 


were  chosen  to  forestall  disproportionate  contributions  from 
large  c.m.  separations  (unlikely  to  be  of  much  importance  in 
the  dense  fluid  regime  sampled  below)  by  setting  <psw 
=  1000  K,  and  the  weighted  root  mean  square  error 
(wRMSE)  was  evaluated  as 

/  Err 

wRMSE=  a/ =4—.  (11) 

v  ^1=1^1 


This  procedure  yielded  the  following  best-fit  parameter  val¬ 
ues 


e  =  34.156  K, 

r0  =  4.037  A,  (12) 

a=  12.29, 

with  a  wRMSE  of  0.189.  At  fixed  parameter  values  <pab  is 
determined  entirely  by  the  site-site  separation  (rab),  but  it  is 
important  to  note  that  the  full  interaction  of  two  diatomics 
((Pij)  is  highly  anisotropic.  The  character  of  the  reference 
potential  defined  by  Eqs.  (6)-(8)  and  (12)  is  illustrated  in 
Fig.  2,  where  we  have  plotted  over  a  range  of  r(/  for  the 
four  fiducial  configuration  types  shown  in  Fig.  1.  These 
should  be  compared  not  only  with  one  another  but  also  with 
a  purely  spherical  potential  (/=0)  for  which  the  four  configu¬ 
ration  types  collapse  into  one. 

In  order  to  estimate  the  magnitude  of  many-body  effects, 
we  compared  the  standard  pair  energy, 

U2=UAB-UA-UB,  (13) 

with  an  “effective”  pair  energy  for  molecules  in  the  liquid, 

U AB  ~  ^b°x  ^bo x-AB  ~  U W)x- A  ~  ^box-B  ’  ( 1  ^) 

where  U A  ( UB)  is  the  total  energy  of  isolated  molecule  A(B), 
UAB  is  that  of  two  molecules  taken  together,  f/box  is  the  en¬ 


ergy  of  a  box  containing  100  molecules  in  the  liquid  state, 
and  Ubm_x  is  that  of  the  same  box  with  the  designated  mol¬ 
ecule  or  pair  removed.  Under  the  assumption  that  many-body 
effects  terminate  at  three-body  interactions,  it  can  be  shown 
(see  Appendix  A)  that  the  sum  over  all  A  and  B  of  U* 

-  ififl  gives  the  many-body  energy  of  the  sample  to  within  a 
constant  factor.  For  a  pair  of  molecules  taken  from  a  liquid- 
state  configuration  at  7  =  3267  K  and  P=  16.8  GPa,  UAB/kB 
=  4519  K  and  U*AB/kB= 3030  K  [it  is  important  to  note  that  in 
evaluating  Eqs.  (13)  and  (14),  all  atoms  removed  were  re¬ 
placed  with  a  ghost  basis  set  just  as  in  standard  counterpoise 
calculations],  giving  a  many-body  stabilization  per  pair  of 
approximately  1500  K,  or  one-third  of  the  isolated  pair  en¬ 
ergy.  Although  a  rough  estimate  only,  this  provides  strong 
prima  facie  evidence  of  an  appreciable  many-body  contribu¬ 
tion  to  the  energy,  and  thus  motivation  for  moving  beyond 
straightforward  use  of  pair  potentials  when  high  accuracy  is 
required. 


III.  NESTED  MARKOV  CHAIN  MONTE  CARLO 
SAMPLING 


The  following  provides  only  a  brief  overview  of 
o-n(MC)2;  for  a  more  complete  treatment,  see  Paper  1.  In 
standard  MC  sampling  with  the  Metropolis  algorithm,  a  trial 
step  drawn  from  a  uniform  distribution  (also  known  as  the 
marginal  distribution)  is  accepted  with  probability33 

/  77/  \ 

a7;-  =  min^  ,lj,  (15) 

where  7 rk  represents  the  relative  weight  assigned  to  configu¬ 
ration  space  point  k.  Trk°^eWk  in  classical  statistics  and  Wk  is 
a  thermodynamic  function  appropriate  to  the  ensemble  being 
sampled.  In  the  isothermal-isobaric  (. NPT )  ensemble  used 
below,  W  is  defined  as 

Wk  =  -p{Uk  +  PVk)  +  N  lnVt,  (16) 


the  final  term  of  which  results  from  the  use  of  scaled 
coordinates.5'34  U,  V,  and  /3  represent  internal  energy,  vol¬ 
ume,  and  inverse  temperature  ( kT)~l ,  respectively.  Trial 
moves  consist  of  single -particle  displacements  and  volume 
adjustments,  and  the  relative  probabilities  of  attempted  move 
types  are  chosen  a  priori  according  to  the  needs  of  the  simu¬ 
lation.  Ensemble  averages  of  some  property  X  can  be  ap¬ 
proximated  by  discrete  quadrature  according  to 


,  ,  fTT(T,V)XdTdV  1  ^ 


(17) 


where  t  represents  the  full  array  of  system  coordinates  and 
the  sum  runs  over  the  Ns  configurations  sampled  by  the  MC 
simulation.  The  variance  in  such  averages  will  decrease  as  Ns 
rises  and,  in  the  limit  of  fully  decorrelated  samples,  the  scal¬ 
ing  relation  assumes  the  well-known  form  crx~N~s112.  Con¬ 
figurations  sampled  consecutively  using  standard  MC  are  re¬ 
lated  by  elementary  moves;  that  is,  they  are  highly 
correlated.  If  the  system  energy  is  expressible  in  a  simple 
analytical  form,  then  the  use  of  highly  correlated  samples 
can  be  compensated  for  by  accumulating  sufficiently  large 
values  of  Ns.  The  difficulty  arises  when  one  requires  both 
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high  accuracy,  such  as  that  offered  by  an  ah  initio  or  DFT 
potential,  and  high  precision,  which  requires  many  decorre- 
lated  sampling  points.  If  each  acceptance  test  (energy  evalu¬ 
ation)  is  costly  in  terms  of  computing  time,  then  being  forced 
to  collect  strongly  correlated  samples  may  preclude  precision 
sampling  on  purely  practical  grounds. 

A  method  for  sampling  an  accurate  and  expensive  poten¬ 
tial  without  requiring  the  strongly  correlated  samples  of  con¬ 
ventional  MC  was  introduced  by  Iftimie  et  al.,9  reformulated 
by  Gelb,10  and  refined  by  the  current  authors.8  The  n(MC)2 
method  concatenates  a  sequence  of  elementary  moves,  each 
accepted  with  Boltzmann  weight  on  the  basis  on  an  approxi¬ 
mate  reference  potential.  At  the  endpoints  of  this  sequence  a 
modified  Metropolis  test  is  used  to  recover  a  Boltzmann  dis¬ 
tribution  of  states  based  on  a  different  (and  hopefully,  more 
accurate)  potential.  We  will  refer  to  the  elementary  steps 
taken  independently  as  reference  steps  (taken  in  the  reference 
system),  and  the  entire  sequence  of  reference  steps  book- 
ended  by  evaluations  of  the  second  potential  as  a  full  step 
(taken  in  the  full  system,  characterized  by  the  full  potential.) 
Consecutive  full  system  energies  are  partially  decorrelated 
by  the  intervening  sequence  of  reference  steps.  If  the  full 
potential  is  more  expensive  than  the  reference  potential,  then 
n(MC)2  lowers  the  total  cost  of  the  calculation  by  limiting 
the  Ns  required  for  a  target  <rx  on  averages  taken  in  the  full 
system.  Indicating  reference  system  quantities  with  super¬ 
scripted  “0,”  the  form  of  the  acceptance  probability  required 
for  Metropolis  sampling  of  the  partially  decorrelated  (full) 
potential  is 

(  nM0)\  ,  , 

^  =  min^,^-pj.  (18) 

To  summarize,  a  sequence  of  trial  steps  taken  in  the  refer¬ 
ence  system  and  accepted  according  to  (15)  is  used  to  build 
a  trial  step  taken  in  the  full  system  that  is  accepted  according 
to  (18).  Note  that  the  standard  ratio  of  Boltzmann  factors  for 
the  initial  and  final  states  in  the  full  system  has  been  “cor¬ 
rected”  by  the  inverse  of  this  ratio  for  the  reference  system. 
Following  our  previous  treatment,  we  express  the  difference 
in  full  and  reference  system  weights  for  state  k  as  SWk 
=  W i  -  W°°  and  the  difference  in  SW  for  consecutive  full 
system  configurations  as  AW=  SWj-SWj.  Equation  (18)  can 
then  be  recast  as 


fl,  AW  3=  0 

a‘J~\em,  AW  <  0. 


(19) 


The  total  speedup  factor  S  gained  by  sampling  the  full 
potential  with  n(MC)2  rather  than  with  traditional  MC  can  be 
given  in  terms  of  the  computing  time  needed  to  reach  a  target 
variance  requiring  Nrl  decorrelated  evaluations  of  the  full  po¬ 
tential.  The  cost  of  collecting  these  directly  through  standard 
Metropolis  sampling  is 


A  -  NdOcottK, 


(20) 


where  X  is  the  computational  expense  of  making  an  elemen¬ 
tary  move  with  the  full  potential.  Ocorr  is  an  effective  corre¬ 
lation  length  (assumed  to  be  roughly  equal  for  reference  and 
full  potentials)  statistically  equivalent  to  a  single  random 


sample  drawn  from  the  distribution.  It  is  only  an  effective 
length  in  the  sense  that  it  incorporates  the  inefficiency  of 
rejected  steps  as  well  as  net  motion  through  configuration 
space,  and  therefore  it  will  vary  with  the  maximum  trial  step 
radius.  For  the  purposes  of  this  discussion,  we  assume  this 
radius  to  be  set  such  that  the  mean  acceptance  probability  of 
elementary  trial  steps  made  using  either  potential  is  the 
same.  The  total  cost  of  collecting  Nd  decorrelated  samples 
using  the  n(MC)2  procedure  is  then 

.,  ,  (oi  NdOcon 

ArfOcorrX  +  X 

A'  = - z - ,  (21) 

a 

where  X(0)  is  the  cost  of  an  elementary  move  made  with  the 
reference  potential,  a  is  the  mean  acceptance  probability  of 
trial  composite  steps,  and  O  is  the  number  of  reference  steps 
comprising  a  single  composite  step.  The  speedup  factor  is  the 
ratio  of  A  and  A', 


aX 

X®+- 

o 


(22) 


In  his  presentation  of  n(MC)2,  Gelb10  made  the  additional 
simplification  that  k/\^-0(N)  when  a  pair  potential  is  used 
as  reference  for  an  /V-body  potential  such  as  DFT.  Because 
the  prefactor  in  such  scaling  relations  can  be  highly  signifi¬ 
cant,  we  include  it  explicitly  and  simplify  Eq.  (22)  further 
such  that 


kN+  O 

where  X/\^-kN.  Although  we  have  not  attempted  to  quan¬ 
tify  k  carefully,  a  conservative  estimate  of  its  value  in  our 
simulations  is  0(  107),  whereas  a  large  value  for  O  would  be 
<9(103).  Therefore  we  ignore  the  first  term  in  the  denomina¬ 
tor  and  obtain  as  an  approximate  speedup  factor 

S  ~  aO.  (24) 

This  is  an  intuitively  satisfying  result,  indicating  that  when 
the  cost  of  evaluating  the  reference  potential  is  negligible 
compared  to  that  of  evaluating  the  full  potential,  the  speedup 
from  using  n(MC)2  over  MC  is  roughly  equal  to  the  effective 
length  (actual  length  tempered  by  average  probability  of  ac¬ 
ceptance)  of  composite  steps  taken  in  the  full  system. 

The  presence  of  a  in  the  denominator  of  Eq.  (21)  sug¬ 
gests  a  simple  prescription  for  raising  the  value  of  S:  In¬ 
crease  the  mean  acceptance  probability  of  trial  composite 
steps.  Although  composite  steps  built  from  a  sequence  of 
elementary  steps  taken  on  the  reference  potential  lower  the 
correlation  between  consecutive  evaluations  of  the  full  po¬ 
tential,  differences  between  the  reference  and  full  system  dis¬ 
tributions  impose  practical  limits  on  the  length  of  these  se¬ 
quences.  That  is,  because  the  distributions  generally  will 
peak  in  different  regions  of  configuration  space,  the  greater 
the  number  of  elementary  steps  taken  on  the  reference  po¬ 
tential  between  acceptance  tests  taken  on  the  full  potential. 
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the  more  the  reference  system  Markov  chain  will  “drift” 
back  to  the  center  of  its  own  distribution  (see  especially  Figs. 
11-14  in  Paper  1)  and  thereby  lower  the  mean  acceptance 
probability  of  trial  composite  steps.  Conversely,  the  greater 
the  difference  in  distribution  of  reference  and  full  potentials, 
the  lower  the  mean  acceptance  probability  of  trial  composite 
steps  built  from  a  fixed  number  of  elementary  steps  taken  on 
the  reference  potential.  It  would  seem,  then,  that  raising  the 
mean  acceptance  probability  of  trial  composite  steps  requires 
shifting  one  or  both  of  the  reference  and  full  distributions 
such  that  their  overlap  is  enhanced. 

As  detailed  in  Paper  1,  the  average  acceptance  probabil¬ 
ity  for  n(MC)2  steps  can  be  expressed  exactly  in  the  limit 
that  a  pair  of  configurations  i  and  j  are  decorrelated.  We 
express  the  mean  acceptance  probability  in  this  limit  as 


A  = 


lim 

o^ocm 


(25) 


An  initial  state  i  will,  by  construction,  possess  relative 
weight  ew‘  if  the  goal  is  to  sample  the  full  distribution.  The 
final  state  j  will,  in  the  OCOTr  limit,  be  drawn  randomly  from 

the  reference  distribution  and  thus  carry  weight  e  /  .  Note 
that  the  acceptance  probability  atj  given  in  Eq.  (18)  also  can 
be  expressed  in  terms  of  SW  and  U/,0)  for  states  i  and  j  so 
that  ajj  averaged  over  the  configuration  and  volume  spaces 
of  all  decorrelated  (ij)  pairs  gives 

-_  ffff  a,iew/+wlj  ]d TjdVjd TjdV , 

ffffeW'+w<J  ] dTtdV jdTjdV j 

_  ffff  a  if  m'(e  wf)+wj>)  d  TjdVjd  TjdVj 

ffffe^fe^  >+wJ  ’jdr/dV/d TjdVj 


The  integrals  in  Eq.  (26)  can  be  estimated  by  sampling  SW  at 
A'rw  configurations  drawn  from  the  reference  distribution;  for 
the  combination  of  potentials  used  here,  Nm=  0(  1 00).  We 
refer  to  this  process  as  the  reweighting  calculation,  and 
estimates  of  A  built  in  this  manner  will  be  denoted  ,4rw.  Note 
that  the  terms  appearing  in  parenthesis  in  Eq.  (26)  are  present 
implicitly  when  i  and  j  are  sampled  from  the  reference  dis¬ 
tribution  so  that  Eq.  (26)  can  be  written  more  succinctly  as 


-  {{a,?*”))  o 

«eW»o 


(27) 


where  the  subscripted  zeros  indicate  that  the  distribution  be¬ 
ing  sampled  is  that  of  the  reference  potential.  Because  the 
SW  appearing  in  Eq.  (26)  depend  on  the  full  system  pressure 
and  temperature  through  Eq.  (16),  a  single  set  of  /Vrw  points 
sampled  from  the  reference  distribution  at  a  given  Pi{))  and 
1<<]>  can  be  used  to  build  a  family  of  Arw  varying  parametri¬ 
cally  as  functions  of  P  and  T.  It  is  appropriate  to  express  Arw 
as  a  function  of  all  four  thermodynamic  variables,  Alw 
= yf/4°), 7<°), /j, y-j,  anc|  ajj  subsequent  references  to  Arw  will 
list  them  in  this  order.  Note  that  Arw  is  an  a  priori  estimate  of 
acceptance  probability  for  composite  steps  taken  with  the 
n(MC)2  procedure  in  that  it  predicts  acceptance  probabilities 
without  reliance  upon  an  actual  n(MC)2  simulation.  If  one 
desires  full  system  thermodynamic  data  at  ( P=P ',  T=T'),  a 


reweighting  calculation  performed  at  those  conditions  can  be 
followed  by  variation  in  P  and  T  so  as  to  maximize  Arw.  We 
express  this  optimization  procedure  as  follows: 

« (max)  .  r  pi  t'  p  T  \ 

-^rw  •> 1  ■> r  opt’ 1  opt/ 

=  max{Arw(F(0),  7<0\x,y):P{0)  =  P',Ti°^  =  T':Pl 

=S  x  P2,T|  =£  y  =£  T2}.  (28) 

The  initial  domain  covered  by  the  full  system  variables  can 
be  chosen  generously,  then  squeezed  incrementally  in  con¬ 
junction  with  finer  meshes;  we  have  found  the  Arw  surface  to 
vary  smoothly  with  (x  ,y)  =  (P  ,T),  as  illustrated  in  Fig.  6  of 
Paper  1.  There  we  obtained  (P=Popt,  7  =  7),pl)  for  full  system 
variables,  then  extrapolated  linearly  back  to  the  original  (P 
=  P',  T=T'),  and  applied  the  same  transformation  to  the  ref¬ 
erence  system  variables  to  yield  (P®  ~  ~  7^) .  In 

this  way,  optimized  n(MC)2  simulations  could  be  carried  out 
at  the  full  system  conditions  desired  but  with  reference  sys¬ 
tem  conditions  designed  to  maximize  the  a  priori  acceptance 
probability  of  composite  steps  taken  in  the  full  system.  The 
approach  taken  here  is  slightly  different  in  that  we  performed 
n(MC)2  sampling  with  the  optimal  set  (P(0,  =  P\  7^0)  =  r', 
P=Popt,  7  =  7  opi)  resulting  from  Eq.  (28).  This  generated  full 
system  thermodynamic  conditions  different  from  the  original 
(P' ,T'),  and  an  expansion  procedure  described  in  the  next 
section  was  used  to  recover  results  at  (P,  T)  on  the  Hugoniot. 

Data  illustrating  the  reweighting,  optimization,  and  opti¬ 
mized  n(MC)2  simulation  results  are  shown  in  Table  III.  The 
SW  distribution  was  sampled  at  /Vrw  configurations  drawn 
from  the  reference  distribution  at  the  (P(0)  =  P\  7i(ll=  7’') 
pairs  shown  in  the  two  leftmost  columns.  a|°* 
=Arw(P' ,  T  ,P' ,  T')  was  built  using  Eq.  (27),  then  maxi¬ 
mized  according  to  Eq.  (28)  to  yield  A|™ax)  at  (l’=  Pop[,  T 
=  7'opt);  the  resulting  Popt  and  Popt  are  shown  in  the  third  and 
fourth  columns  of  Table  III.  One  measure  of  the  procedure’s 
effectiveness  is  the  ratio  of  the  optimized  to  original  a  priori 
acceptance  probabilities,  A(max)/A<0>;  as  shown  in  the  right 
side  of  the  table,  the  predicted  acceptance  probabilities  rose 
by  factors  of  roughly  1.5-25.  Note  that  minor  absolute  varia¬ 
tions  in  A^J*  (a  small  number)  can  make  large  differences  in 
the  ratio.  After  maximizing  Arw  through  variation  in  P  and  T 
as  in  Eq.  (28),  optimized  n(MC)2  simulations  were  per¬ 
formed  using  the  (P' ,  P' ,Popt,  Popt)  combinations  that  re¬ 
sulted.  The  number  of  elementary  steps  taken  in  the  refer¬ 
ence  system  used  to  build  a  composite  step  made  in  the  full 
system  was  (9=150,  whereas  (9corr~1000  for  the  reference 
and  full  potentials  chosen  in  Sec.  II  (cf.  Fig.  3  in  Paper  1). 
The  average  speedup  factor  S  achieved  in  sampling  the  full 
potential  with  o-n(MC)2  rather  than  traditional  MC  was  just 
over  50. 

No  attempt  was  made  in  the  present  work  to  maximize 
(9,  the  number  of  elementary  steps  in  the  reference  system 
used  to  build  a  trial  composite  step  in  the  full  system.  This 
potential  loss  of  efficiency  is  compensated  for  partly  by  a 
higher  mean  acceptance  probability  at  lower  (9  values,  but 
the  ideal  situation  would  be  one  in  which  optimization  were 
performed  directly  on  S  rather  than  on  Arw.  This,  in  turn, 
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TABLE  III.  Summary  of  the  reweighting  and  optimization  procedures  described  in  the  text.  The  SW  distribution 
was  sampled  at  configurations  defined  by  the  thermodynamic  conditions  (/Jl0'  =  /-*' ,  7101- T')-  The  a  priori 
acceptance  probability  A® =Am(P' ,  T' ,P' .T')  was  then  calculated  from  Eq.  (27),  followed  by  optimization 
according  to  Eq.  (28)  to  yield  =Arw(P'  ,T'  ,Popt,  Topt).  The  effectiveness  of  the  optimization  is  reflected  in 
the  ratio  A)™2*'/ A®.  The  speedup  factor  S —  aO  is  given  for  optimized  n(MC)2  simulations  in  which  0  =  150 
and  Ocorr»  1000. 


p(0)  =  pt 

jio)=pi 

P 

r  opt 

T’opt 

Wrw 

A® 

(%) 

A”“ 

(%) 

A  max  j  a  (0) 
rw  '  rw 

S 

6.30 

728 

4.84 

615 

563 

0.49 

4.14 

8.45 

41.6 

9.70 

1424 

7.48 

1200 

394 

0.49 

8.84 

18.04 

55.3 

11.90 

1893 

9.22 

1606 

478 

0.82 

3.73 

4.55 

60.4 

14.90 

2540 

11.94 

2206 

430 

0.59 

6.78 

11.49 

55.3 

15.80 

2730 

12.52 

2355 

484 

0.34 

7.96 

23.41 

52.9 

19.20 

3490 

16.79 

3267 

532 

2.94 

9.59 

3.26 

64.7 

20.50 

3795 

16.60 

3270 

444 

0.63 

13.01 

20.65 

53.3 

20.60 

3812 

16.93 

3417 

386 

2.25 

8.62 

3.83 

56.7 

23.00 

4375 

19.98 

4079 

380 

1.74 

7.32 

4.21 

27.7 

23.20 

4410 

19.63 

3886 

444 

2.49 

5.89 

2.37 

76.0 

23.30 

4434 

17.58 

3361 

443 

1.84 

3.09 

1.68 

34.5 

23.50 

4475 

21.27 

4308 

441 

3.72 

5.97 

1.60 

55.2 

24.70 

4765 

19.58 

3854 

525 

0.42 

1.50 

3.57 

60.5 

25.20 

4877 

20.75 

4362 

439 

2.30 

3.65 

1.59 

56.5 

27.40 

5421 

22.01 

4769 

426 

0.76 

3.91 

5.14 

60.7 

28.90 

5785 

21.29 

4765 

420 

1.17 

3.20 

2.74 

40.7 

30.00 

6062 

23.06 

5156 

492 

0.67 

2.91 

4.34 

41.2 

40.50 

8847 

34.90 

8308 

429 

5.48 

7.86 

1.43 

40.0 

would  require  an  analytical  expression  for  the  acceptance 
probability  in  Eq.  (27)  in  terms  of  O,  as  well  as  more  exten¬ 
sive  reweighting  calculations  over  variable  O.  Empirical  ex¬ 
amination  of  the  variation  in  a  with  O  can  be  found  in  Paper 
1. 

In  order  to  limit  the  computational  overhead  of  the  re¬ 
weighting  and  optimization  procedures,  it  is  desirable  that 
converge  for  Nrv/  as  small  as  possible;  otherwise,  effi¬ 
ciency  gained  in  using  the  optimized  form  of  n(MC)2  will  be 
lost  in  performing  the  optimization  procedure  itself.  As  dis¬ 
cussed  in  Paper  1  (and  particularly  in  connection  with  Figs.  7 
and  8  there),  Ap™x)  tends  to  converge  faster  than  (Popt,  Ta pt) 
due  to  its  exhibiting  a  broad  flat  peak  as  a  function  of  P  and 
T.  Figure  3  illustrates  the  convergence  behavior  of  A[™dxl  for 
five  pairs  of  (P=P\  T=T')  spanning  the  full  range  of  con¬ 
ditions  included  in  Table  IV.  All  five  states  exhibited  an  ini¬ 
tial  decline  in  A(max)  as  reweighting  samples  were  added  but 
stabilized  roughly  in  the  neighborhood  of  Arw=  100-500. 
While  convergence  was  incomplete  in  most  cases  (particu¬ 
larly  for  the  more  extreme  thermodynamic  states),  it  still  was 
adequate  to  ensure  large  efficiency  increases  relative  to  stan¬ 
dard  n(MC)2  or  conventional  MC  sampling  on  the  basis  of 
the  full  potential. 

IV.  THE  NITROGEN  HUGONIOT 

The  Hugoniot  locus  is  comprised  of  thermodynamic 
states  accessible  by  shock  loading  of  an  initial  state.  In  this 
section  we  describe  application  of  the  o-n(MC)2  method  to 
construction  of  the  Hugoniot  locus  for  N2  fluid  and  compare 
our  predictions  to  experimental  data.  Prior  to  doing  so,  how¬ 


ever,  we  briefly  summarize  the  continuum  theory  of  shock 
waves  as  well  as  the  statistical  approximations  underlying 
atomistic  treatment. 

Consider  a  sample  of  material  having  specific  volume  V0 
(volume  per  unit  mass=l/p0)  and  specific  internal  energy 
E0,  struck  by  a  piston  moving  at  speed  up.  If  cq  is  the  sound 
speed  in  the  material  at  zero  pressure  and  up  <§  c0,  the  impact 
will  create  an  acoustic  wave  moving  through  the  sample  at 
speed  c0;  pistons  having  speeds  of  up~0.0lc0  or  greater  will 
generate  shock  waves  characterized  by  discontinuous  density 
gradients  at  the  wave  front.  Mass,  momentum,  and  energy 
conservation  in  samples  initially  at  rest  at  zero  pressure 


FIG.  3.  Convergence  behavior  of  ATO  with  respect  to  the  number  of  re¬ 
weighted  sampling  points  A^w,  for  five  combinations  of  initial  ( P=P' ,  T 
=  T).  All  combinations  show  an  initial  decrease  in  acceptance  probability, 
but  level  out  at  roughly  Afrw=500.  Arw  is  an  a  priori  estimate  of  a  in  the 
decorrelated  limit  (0=1000),  and  therefore  much  lower  than  the  a  com¬ 
puted  a  posteriori  from  n(MC)2  simulation  at  0=150. 
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TABLE  IV.  Summary  of  simulation  data  used  in  Figs.  4-8. 


p 

(GPa) 

T 

(K) 

V 

(cm3/  g) 

cr(V) 

(cm3/g) 

U 

(kJ/g) 

<r(U) 

(kJ/g) 

E 

(kJ/g) 

AEC 

(kJ/g) 

AT 

(K) 

A  E 
(kJ/g) 

AV 

(cm3/ g) 

4.84 

615 

0.7126 

0.0090 

0.714 

0.029 

1.174 

0.070 

-31 

-0.045 

-0.0102 

7.48 

1  200 

0.6678 

0.0059 

1.107 

0.060 

2.062 

0.121 

-75 

-0.115 

-0.0015 

8.00 

1  172 

0.6590 

0.0074 

1.178 

0.052 

2.107 

-0.012 

9 

0.010 

0.0006 

9.22 

1  606 

0.6535 

0.0068 

1.431 

0.052 

2.762 

0.276 

-224 

-0.236 

-0.0085 

10.00 

1  557 

0.6299 

0.0080 

1.590 

0.061 

2.875 

0.050 

-39 

-0.041 

-0.0017 

11.94 

2  206 

0.6195 

0.0054 

1.819 

0.056 

3.733 

0.273 

-210 

-0.236 

-0.0061 

12.00 

1  957 

0.6112 

0.0070 

1.873 

0.071 

3.543 

0.015 

-13 

-0.014 

-0.0003 

12.52 

2  355 

0.6228 

0.0088 

1.843 

0.059 

3.904 

0.291 

-207 

-0.220 

-0.0114 

15.00 

2  585 

0.5863 

0.0068 

2.337 

0.091 

4.628 

0.001 

-1 

-0.001 

0.0000 

16.60 

3  270 

0.5874 

0.0092 

2.594 

0.094 

5.575 

0.451 

-350 

-0.348 

-0.0123 

16.79 

3  267 

0.5872 

0.0074 

2.530 

0.079 

5.508 

0.322 

-253 

-0.285 

-0.0044 

16.93 

3  417 

0.5760 

0.0060 

2.737 

0.097 

5.867 

0.543 

-464 

-0.536 

-0.0008 

17.58 

3  361 

0.5686 

0.0054 

2.823 

0.148 

5.896 

0.297 

-200 

-0.296 

-0.0001 

19.58 

3  854 

0.5682 

0.0079 

3.096 

0.153 

6.671 

0.417 

-310 

-0.357 

-0.0061 

19.63 

3  886 

0.5594 

0.0066 

3.074 

0.074 

6.681 

0.325 

-250 

-0.282 

-0.0044 

19.98 

4  079 

0.5508 

0.0103 

3.053 

0.088 

6.858 

0.299 

-158 

-0.218 

-0.0081 

20.00 

3  704 

0.5514 

0.0070 

3.137 

0.097 

6.559 

0.001 

0 

-0.001 

0.0000 

20.75 

4  362 

0.5544 

0.0063 

3.307 

0.106 

7.401 

0.623 

-425 

-0.544 

-0.0076 

21.27 

4  308 

0.5583 

0.0070 

3.172 

0.105 

7.211 

0.302 

-229 

-0.266 

-0.0033 

21.29 

4  765 

0.5599 

0.0064 

3.449 

0.104 

7.957 

1.058 

-769 

-0.940 

-0.0111 

22.01 

4  769 

0.5482 

0.0102 

3.439 

0.121 

7.951 

0.686 

-549 

-0.532 

-0.0139 

23.06 

5  156 

0.5404 

0.0087 

3.722 

0.178 

8.631 

0.924 

-728 

-0.852 

-0.0062 

25.00 

4915 

0.5286 

0.0093 

3.935 

0.245 

8.597 

0.083 

-61 

-0.079 

-0.0003 

30.00 

6218 

0.5123 

0.0044 

4.426 

0.108 

10.429 

-0.056 

49 

0.054 

0.0002 

34.90 

8  308 

0.5054 

0.0054 

5.319 

0.183 

13.483 

1.144 

-951 

-1.089 

-0.0031 

35.00 

7  612 

0.4999 

0.0079 

5.194 

0.156 

12.638 

0.166 

-141 

-0.144 

-0.0013 

40.00 

9  096 

0.4775 

0.0079 

6.179 

0.533 

15.158 

0.441 

-290 

-0.466 

0.0012 

45.00 

10  675 

0.4752 

0.0032 

6.387 

0.132 

17.002 

0.377 

-324 

-0.356 

-0.0009 

( Pq~0 )  yield  for  specific  volume  ( VH),  density  (pH=[IVH), 
pressure  ( PH ),  and  specific  internal  energy  (EH)  behind  the 
shock  front36 

PH=p0usup,  (29) 


Ph=Po 


(30) 


Eh  =  E0+1-Ph(V0-Vh)  =  E0+^Ph(^-^),  (31) 

2  2  \  pHp0  J 

where  us  is  the  shock  front  velocity,  up  is  the  particle  velocity 
(equal  to  the  piston  velocity),  and  final  values  have  been 
designated  by  a  subscripted  H  for  Hugoniot.  These  expres¬ 
sions  are  known  collectively  as  the  Rankine-Hugoniot  rela¬ 
tions  or  jump  conditions.  Equations  (29)— (3 1)  identify  the 
locus  of  states  thermodynamically  accessible  to  the  sample 
upon  shock  compression  of  a  given  initial  state,  and  in  this 
they  embed  considerable  information  regarding  the  material 
equation  of  state  (EOS). 

Two  important  variables  in  Eq.  (31)  are  E,  the  specific 
internal  energy  or  energy  per  unit  mass,  and  V,  specific  vol¬ 
ume  or  volume  per  unit  mass.  Thermodynamic  contributions 
and  simulation  averages  are  discussed  most  naturally  in 
terms  of  total  energy  and  volume  of  an  A'-particle  system  or 
even  energy  and  volume  per  particle;  in  order  to  maximize 
transparency,  we  will  employ  the  same  variable  notation  in 


all  these  cases.  That  is,  E  will  designate  specific  internal 
energy  as  well  as  energy  per  particle  and  energy  of  an 
/V-particle  system.  Likewise,  V  will  designate  related  volume 
variables.  It  is  straightforward  to  distinguish  among  usages 
by  context  or  units.  Note  that  the  conversion  factor  MINA  is 
g/particle  where  M  is  the  molecular  weight  ( M 
=  28.0134  g/mol  for  N2)  and  NA  is  Avogadro’s  number 
(6.0221  367  X  1023  particles /mol).  Multiplying  both  sides  of 
Eq.  (31)  by  MINA  changes  the  E  and  V  variables  from 
energy/g  and  volume/g  to  energy/particle  and  volume/ 
particle,  leaving  the  equation  unchanged  except  for  the 
meaning  of  E  and  V. 

E  readily  separates  into  three  distinct  components  for  a 
/V-particle  system, 

E  =  Eid  +  Ev  +  Eex,  (32) 

the  first  two  of  which  can  be  evaluated  using  standard 
forms37 

EiA=\NkBT  (33) 


is  the  ideal  translational  and  rotational  contribution  of  a  lin¬ 
ear  rigid  rotor,  and 


NkB®0 

e&»,T-  1 


(34) 


is  the  energy  of  a  quantum  harmonic  oscillator.  The  tempera¬ 
ture  range  sampled  was  too  low  for  electronic  excitation  to 
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play  a  significant  role,38  and  the  vibrational  temperature  ©„ 
=  ha>/kB  was  evaluated  using  the  experimentally  determined 
nitrogen  frequency  at  zero  pressure.  9  Vibrational  zero  point 
energy  was  included  separately  in  the  reference  heat  of  for¬ 
mation  and  specific  internal  energy  E0,  which  were  deter¬ 
mined  from  measured  thermodynamic  properties  of  liquid 
N2.  The  zero  of  energy  is  that  for  isolated  N2  molecules  at 
T=0  K.  The  excess  contribution  Eex  is  defined  by 

Eex  =  (U),  (35) 

where  ( U )  is  the  average  configurational  energy  calculated 
from  o-n(MC)2  simulation. 

According  to  Eq.  (31),  Eh  corresponds  to  values  of  E 
which  are  zeros  of  the  function 


f{E)  =  E-E0-\p{V 0-VH), 


(36) 


where  VH  is  chosen  as  an  independent  variable,  P(VH,E)  is 
determined  by  the  EOS,  and  EH(VH )  is  the  dependent  vari¬ 
able  expressed  as  a  function  of  V H.  Solutions  to  Eq.  (36) 
define  a  locus  of  pressure  and  temperature  combinations  (P 
=  PH ,  T=Th)  comprising  the  Hugoniot.  Maximization  of  ,4 rw 
with  Eq.  (28)  generated  a  set  of  ( P=Popt ,  7  =  /„pt)  at  which 
o-n(MC)2  simulations  were  performed;  this  set  differed  from 
the  original,  desired  (P=P',  7  =  T').  The  original  combina¬ 
tions  were  chosen  on  the  basis  of  previously  reported  re¬ 
weighting  calculations,40’41  as  states  near  the  Hugoniot  locus. 
In  order  to  extract  thermodynamic  information  at  Hugoniot 
states  from  o-n(MC)2  results  collected  at  ( P ,  7)  pairs  not  on 
the  Hugoniot,  a  mapping  between  the  two  sets  of  conditions 
must  be  provided.  For  (P,T)  sufficiently  close  to  ( PH ,  Th), 
volume  and  energy  can  be  expanded  to  first  order  about  their 
simulation  values. 


Vh-V=AV~ 


dV 

dT 


(Th-T)  + 


dV 

dP 


e pH-p ) 


dV 

dT 


•  A  T  + 

p 


dV 

dP 


•  A  P, 

T 


(37) 


EH-  E  =  AE  ~ 


dE 

dT 


(T„-T)  + 


dE 

dP 


( PH-P ) 


dE 

dT 


■  AT  + 

p 


dE 

dP 


■  A  P. 

T 


(38) 


Here,  P  and  T  are  fixed  parameters  in  the  NPT  ensemble, 
while  U  and  V  represent  the  average  energy  and  volume 
(U ={U) P  T,  V=(V)P  T)  calculated  from  o-n(MC)2  simulation 
at  a  prescribed  pressure  and  temperature;  E  is  related  to  (U) 
by  Eqs.  (32)-(35).  If  P=Ph  by  construction,  the  second  of 
the  two  terms  in  each  of  the  above  expansions  vanishes;  the 
goal  then  is  to  recover  the  Hugoniot  temperature  Th  corre¬ 
sponding  to  a  given  PH.  Toward  this  end,  Eq.  (31)  can  be 
restated  in  terms  of  differences  between  simulation  results 
and  Hugoniot  values, 

E+AE  =  E0+\ph{V0-V-AV),  (39) 

into  which  substitution  of  Eqs.  (37)  and  (38)  produces 
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FIG.  4.  Partial  derivatives  of  volume  with  respect  to  temperature  at  fixed 
pressure  (left  panel),  as  calculated  from  results  of  o-n(MC)2  simulation  per¬ 
formed  at  the  (P=Popt,  T=Topt)  combinations  shown  in  Table  IV.  The  right 
panel  shows  the  same  data,  but  following  multiplication  by  pressure;  as 
described  in  the  text,  these  pressures  are  taken  to  correspond  to  Hugoniot 
values  ( P=PH )  by  construction.  The  partial  derivative  is  evaluated  using  Eq. 
(44),  and  the  dashed  lines  represent  quadratic  (left)  and  linear  (right)  fits  to 
the  data.  Note  the  difference  in  ordinate  scales. 


E  + 


dE 

dT 


•A  T=E0+l-PH[v0-V- 


(40) 


Using  Eq.  (31)  to  define  the  difference  between  the  com¬ 
puted  energy  E  and  EH,  then  substituting  in  Eq.  (40)  gives 


A  EC  =  E-E0- 


1 

2 


PH(V o-V) 


which  yields  upon  slight  rearrangement 


AT  = 


-A  Ec 


/  dE 

1 

+  ~PH 

dV 

) 

\  dT 

P  L 

dT 

P* 

(41) 


(42) 


the  desired  relation.  The  ideal  and  quantum  harmonic  oscil¬ 
lator  contributions  to  dEI dT  are  trivial  to  evaluate  using  Eqs. 
(33)  and  (34),  and  the  rest  of  the  denominator  in  Eq.  (42)  can 
easily  be  recast  (see  Appendix  B)  in  terms  of  simulation 
parameters  and  results 


d(U) 

dT 


=  £[P«UV)  -  (U)(V))  +  ( U 2>  -  (U)2l 


(43) 


djV) 

dT 


p 


~,[p((uv)  -  (U)(V))  +  (V2)  -  (V)2], 


(44) 


where  ensemble  averages  once  again  have  been  indicated 
explicitly.  PH  and  Th  can  now  be  substituted  into  Eqs.  (37) 
and  (38)  to  produce  Eh  and  VH-  Upon  combination  of  Eqs. 
(29)  and  (30)  to  generate  us  and  up,  the  Hugoniot  locus  is 
completely  specified. 

Figures  4  and  5  illustrate  in  more  detail  the  manner  in 
which  Eq.  (42)  was  used  to  recover  temperature  values  on 
the  Hugoniot.  A Ec  was  evaluated  using  Eqs.  (33)-(35),  with 
(U)  given  by  the  average  DFT  energy  predicted  from 
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FIG.  5.  Evaluation  of  the  denominator  appearing  in  Eq.  (42)  as  calculated 
from  results  of  o-n(MC)2  simulation  performed  at  the  (P=Popt,  7’=ropt) 
combinations  shown  in  Table  IV.  As  described  in  the  text,  these  pressures 
are  taken  to  correspond  to  Hugoniot  values  ( P=PH )  by  construction  and  the 
corresponding  Hugoniot  values  for  temperature  are  recovered  by  use  of  Eq. 
(42).  The  dashed  line  is  a  linear  fit  to  the  data. 

o-n(MC)2  simulation  at  the  pressure  and  temperature  combi¬ 
nations  shown  in  Table  III.  The  pressures  at  which  simula¬ 
tions  were  performed  were  also  taken  to  fall  on  the  Hugoniot 
(vide  supra),  and  the  derivatives  appearing  in  the  denomina¬ 
tor  of  Eq.  (42)  were  evaluated  using  Eqs.  (43)  and  (44)  and 
the  equilibrium  averages  computed  from  simulation.  The  de¬ 
rivative  of  volume  with  respect  to  temperature  at  each  ( P 
=  /Jop,  =  PH,  /  =  /opi)  combination  is  shown  in  Fig.  4  (left 
panel),  as  well  these  derivatives  multiplied  by  their  corre¬ 
sponding  pressures  (right  panel).  The  data  in  the  left  and 
right  panels  are  fit  to  linear  and  quadratic  functions  (dashed 
lines),  respectively.  Figure  5  illustrates  the  entire  denomina¬ 
tor  of  Eq.  (42)  built  from  the  component  pieces  appearing  in 
Fig.  4. 

Recovery  of  Th  permitted  determination  of  VH  through 
application  of  Eq.  (37),  enabling  construction  of  the  N2 
Hugoniot  locus  in  Fig.  6  (P-V  plane)  and  Fig.  7  (us-up 
plane).  The  nitrogen  Hugoniot  has  been  measured4"-44  or 
calculated38,4041  on  a  number  of  occasions  previously,  but 
only  in  a  few  instances  have  calculations  using  a  self- 
consistent  potential  such  as  DFT40’45,46  been  reported.  Open 
circles  in  Fig.  6  represent  o-n(MC)2  results  calculated  at 
(Popt, 70pi)  and  then  shifted  to  the  Hugoniot  locus  ( PH,TH )  in 
the  manner  described  above.  The  solid  line  is  a  fit  of  those 
results  to  an  exponential  function.  Although  components  of 
the  shifts  (illustrated  in  Figs.  4  and  5)  possess  relatively  large 
uncertainties,  the  net  uncertainties  are  comparable  to  or 
smaller  than  the  statistical  uncertainties  in  the  simulated 
states  themselves;  error  estimates  are  given  in  Table  IV  but 
omitted  from  the  figures  for  clarity.  The  28  MC  data  points 
shown  explicitly  in  Fig.  6  have  been  replaced  with  a  qua¬ 
dratic  fit  in  Fig.  7  for  ease  of  viewing,  but  otherwise  symbols 
and  shading  have  been  preserved  across  the  two  figures  for 
each  data  set  depicted.  Agreement  with  experimental  data  in 
the  molecular  regime — which  we  take  to  be  40-45  GPa  or 
less — taken  from  Refs.  42  (squares),  43  (filled  circles),  and 
44  (triangles)  is  excellent.  [See  the  following  paragraph  for 
discussion  of  the  higher  pressure  region  for  which  the 


V  /(cm3/g) 

FIG.  6.  N2  Hugoniot  in  the  P-V  plane.  Experimental  results  (filled  symbols) 
are  taken  from  Refs.  42-44;  theoretical  points  were  calculated  using 
o-N(MC)2  (open  circles,  present  work)  and  AIMD  (cross.  Ref.  45).  MC 
simulations  were  restricted  to  the  molecular  regime  (up  to  P=40-45  GPa), 
and  have  been  shifted  according  to  Eq.  (42).  The  inset  displays  the  high- 
pressure  region  that  includes  ionization;  an  AIMD  result  in  the  transition 
region  has  been  circled  in  both  the  inset  and  main  figure  in  order  to  facilitate 
visual  passage  between  the  two.  The  solid  line  is  an  exponential  fit  through 
the  28  o-N(MC)2  data  points. 

o-n(MC)2  methodology  was  not  applied.]  In  addition  to  the 
experimental  results  we  have  included  those  of  a  previous 
calculation43  made  using  ab  initio  molecular  dynamics 
(AIMD)  in  the  NVT  ensemble  with  a  PW91  exchange- 
correlation  functional  and  Vanderbilt  ultrasoft 
pseudopotentials47  (crosses).  This  particular  implementation 
was  designed  specifically  to  capture  electronic  and  dissocia¬ 
tion  effects  in  the  high  temperature  and  pressure  region; 
however,  the  small  system  size  (32  atoms  =  16  molecules) 
and  the  use  of  pseudopotentials  introduce  uncertainties  that 
are  difficult  to  quantify.  Only  four  of  the  calculated  points 
fell  in  the  molecular  regime,  but  these  were  in  qualitative 
agreement  with  experiment.  The  shock  and  particle  veloci- 


FIG.  7.  N2  Hugoniot  in  the  us-up  plane.  Experimental  results  (filled  sym¬ 
bols)  are  taken  from  Refs.  42-44;  theoretical  results  were  calculated  using 
o-N(MC)2  (solid  line,  present  work)  and  AIMD  (cross  Ref.  45).  MC  simu¬ 
lations  were  restricted  to  the  molecular  regime  (up  to  P=40  45  GPa),  and 
have  been  shifted  according  to  Eq.  (42).  Here,  for  ease  of  viewing,  the  28 
individual  o-N(MC)2  points  results  shown  explicitly  in  Fig.  6  were  replaced 
by  a  quadratic  fit  after  transformation  between  the  P-V  and  us-up  planes. 
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ties,  us  and  up,  in  Fig.  7  are  based  on  the  same  fundamental 
information  as  are  P  and  V  in  Fig.  6  {PH  and  VH  are  derived 
from  us,  up,  and  the  Hugoniot  jump  conditions,  Eqs. 
(29)— (3 1),  although  the  relationship  in  the  us-up  plane  is 
nearly  linear.  Again,  however,  the  trend  is  that  the  present 
calculations  are  in  excellent  agreement  with  the  experimental 
data  below  up~ 6  krn/s. 

Thus  far  we  have  focused  exclusively  on  pressures  of 
45  GPa  and  below,  but  Figs.  6  and  7  also  include  higher 
pressure  regions  where  dissociation,  ionization,  and  elec¬ 
tronically  excited  states  become  important.  Trends  in  the  ex¬ 
perimental  data  of  Nellis  et  a/.4,  and  AIMD  simulation  re¬ 
sults  of  Kress  et  a/.43  are  in  good  agreement  generally,  and 
both  clearly  reflect  an  abrupt  change  between  P=40  and 
47  GPa;  here  we  focus  specifically  on  the  transition  region  in 
which  dissociation  emerges  as  an  important  factor.  At  P 
=  31  GPa,  the  specific  volume  predicted  by  ab  initio  molecu¬ 
lar  dynamics  AIMD  (with  7  ±  3%  dissociation)  is  near  the 
average  of  the  measured  data  scatter  (Fig.  6);  at  39  GPa 
(16  ±9%  dissociation  by  interpolation  of  Table  I  in  Ref.  45) 
it  is  significantly  less  than  that  predicted  by  experiment. 
These  facts  suggest  that  the  calculated  dissociation  is  over¬ 
estimated  in  the  transition  region.  By  way  of  contrast,  the 
present  results  fall  well  within  the  experimental  scatter  for 
pressures  up  to  40  GPa.  If  dissociation  were  important  in  this 
region,  we  would  expect  the  o-N(MC)2  results  to  begin  to 
deviate  from  experiment  at  pressures  lower  than  those  for 
which  they  actually  do.  Taken  collectively,  the  measured  and 
calculated  data  suggest  that  dissociation  is  not  significant  to 
the  Hugoniot  below  the  transition  region  found  at  P 
—  40-45  GPa.  This  transition  region  could  reflect  to  a  two- 
wave  structure,  in  which  case  shock  velocity  us  is  established 
by  a  leading  wave  whereas  the  particle  velocity  up  is  deter¬ 
mined  on  the  basis  of  either  thermodynamics  or  kinetics  that 
leads  to  a  secondary  wave  propagating  in  the  material.  Ex¬ 
amples  where  this  could  occur  include  a  first  order  phase 
transition  with  a  large  volume  change  or  relaxation  from  the 
highest  metastable  state  that  can  be  maintained  before  the 
kinetics  allow  a  return  to  equilibrium  on  a  time  scale  that  is 
short  compared  to  the  time  resolution  of  the  shock  experi¬ 
ment.  If  a  single  wave  is  assumed  in  such  an  instance,  then 
the  predicted  P  and  V  will  be  incorrect.  The  hallmark  of 
two-wave  structure  is  constant  us  across  some  interval  of  up\ 
but,  unfortunately,  there  is  only  one  data  point  in  the  transi¬ 
tion  region  shown.  For  this  reason,  it  is  unclear  whether  the 
data  reflect  a  bona  fide  two-wave  structure  or  merely  a  region 
of  rapid  transition. 

Note  that  the  present  calculations  assume  the  validity  of 
the  rigid  rotor  and  quantum  harmonic  oscillator  approxima¬ 
tions.  Raman  line  shifts  in  liquid  nitrogen  singly  shocked  to 
P=  10-20  GPa  and  doubly  shocked  to  P=  1 5-40  GPa  have 
been  measured44  and  calculated  using  a  first-principles 
approach.49  Both  fundamental  (0— >1)  and  hot  band  (1— >2, 
2— >3,  and  3  —  4)  transitions  indicated  that  vibrational  ener¬ 
gies  increased  —1%  in  the  singly  shocked  state  and  —2%  in 
the  doubly  shocked  state.  There  was  no  indication  of  vibra¬ 
tional  potential  softening  even  at  the  highest  pressure  mea¬ 
sured;  the  data  imply  rather  that  vibrations  are  very  weakly 
coupled  to  other  modes,  in  conformity  with  the  assumptions 


FIG.  8.  Temperatures  measured  (filled  symbols,  Refs.  43,  44,  and  48)  or 
calculated  (Ref.  45  and  present  work)  on  the  N2  Hugoniot. 

made  here.  We  have  treated  shocked  nitrogen  in  the  molecu¬ 
lar  phase  as  undergoing  negligible  electronic  excitation  as 
well,  given  that  the  lowest  electronic  excited  state  in  N2  has 
a  relative  energy  of  roughly  6  eV.  This  energy,  and  the 
ground  state  dissociation  energy  of  approximately  10  eV,  are 
to  be  compared  with  a  temperature  equivalent  of  1  eV  or  less 
on  the  principal  Hugoniot  up  to  P  —  45  GPa. 

Several  techniques  for  temperature  measurement  on  the 
Hugoniot  have  been  reported  previously,  each  of  which  po¬ 
tentially  introduces  systematic  errors  of  unknown  magnitude 
due  to  assumptions  made  in  the  analysis.  Figure  8  shows 
measured  data  compared  with  the  present  calculations  and 
those  of  Kress  et  a/.43  The  temperature  data  of  Nellis  et  al ,43 
increase  along  a  trend  line  that  coincides  well  with  the  MC 
results  up  to  P= 36  GPa,  as  would  be  expected  under  condi¬ 
tions  of  negligible  dissociation.  The  next  highest  temperature 
measurement,  obtained  at  P=47  GPa  is  significantly  cooler 
than  this  trend  predicts,  as  would  be  expected  should  onset  of 
significant  dissociation  occur  for  conditions  between  the  two 
measured  data  points.  A  new  trend  is  established  at  higher 
pressure,  corresponding  to  a  partially  dissociated  phase  in 
which  some  fraction  of  kinetic  energy  has  been  expended  to 
break  the  nitrogen  triple  bond,  resulting  in  lower  tempera¬ 
tures  than  would  be  found  for  molecular  nitrogen.  The  inter¬ 
pretation  of  the  temperature  data  is  the  same  as  that  of  the 
P-V  and  us-up  results:  There  is  negligible  dissociation  at 
pressures  below  40-45  GPa,  but  a  significant  amount  at 
pressures  higher  than  this. 

V.  SUMMARY  AND  CONCLUSIONS 

We  have  used  the  o-n(MC)2  algorithm  with  a  DFT  po¬ 
tential  to  generate  the  Hugoniot  locus  of  N2  fluid  for  pres¬ 
sures  and  temperatures  relevant  to  HE  detonation.  After  col¬ 
lecting  several  hundred  configurations  at  which  to  evaluate 
SW  for  a  number  of  ( P ,  7  j  combinations,  the  a  priori  accep¬ 
tance  probability  defined  by  Eq.  (27)  was  maximized  as  a 
function  of  P  and  T  to  yield  optimal  thermodynamic  condi¬ 
tions  at  which  to  perform  n(MC)2  sampling.  The  o-n(MC)2 
simulations  were  then  carried  out  at  temperature  and  pressure 
combinations  not  necessarily  on  the  Hugoniot  locus,  but 
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from  which  Hugoniot  states  could  be  recovered  by  use  of  the 
expansion  procedure  defined  in  Eqs.  (37)  and  (38)  and  cul¬ 
minating  in  Eqs.  (42)-(44).  The  method  yields  both  high 
accuracy  and  high  precision.  It  is  worth  noting  that  in  con¬ 
trast  to  previous  calculations45  the  current  approach  includes 
all  electrons  (no  pseudopotential  approximation)  and  readily 
accommodates  a  system  more  than  six  times  larger  than  the 
earlier  study,  thereby  reducing  the  likelihood  of  finite-size 
effects  in  the  results. 

An  analytical  reference  potential  was  used  in  this  work, 
but  o-n(MC)2  should  work  equally  well  with  two  different 
ab  initio ,  DFT,  or  semiempirical  potentials.  The  reference 
potential  could  rely  on  a  more  approximate  algorithm  or  a 
smaller  basis  set,  so  long  as  (1)  its  correspondence  to  the  full 
potential  still  permits  a  reasonable  offset  O  between  full  en¬ 
ergy  evaluations  and  (2)  it  can  still  be  evaluated  rapidly.  The 
full  potential  would  then  consist  of  a  fully  converged  calcu¬ 
lation  at  a  higher  level  of  theory. 

The  o-n(MC)2  method  provides  a  versatile  framework  in 
which  to  carry  out  precision  MC  sampling  of  accurate  poten¬ 
tials,  and  it  would  be  worthwhile  to  make  detailed  compari¬ 
son  of  its  performance  to  that  of  AIMD.  °  Straightforward 
formulation  of  efficiency  metrics  such  as  Eq.  (24)  is  impos¬ 
sible  without  determination  of  a  quantitative  relation  be¬ 
tween  the  average  single-particle  displacement  in  accepted 
MC  steps  and  in  AIMD  time  steps.  This  is  a  nontrivial  task, 
and  perhaps  it  would  be  simpler  to  calculate  the  number  of 
AIMD  steps  required  to  yield  decorrelated  configurations, 
then  compare  this  with  the  number  of  composite  o-n(MC)2 
steps  needed  for  the  same  purpose;  even  given  these  two 
quantities,  still  it  is  not  at  all  clear  that  their  ratio  would  not 
vary  strongly  from  one  system  to  the  next.  Aside  from  effi¬ 
ciency  considerations,  MC  methods  accommodate  general¬ 
ized  ensemble  sampling  much  more  naturally  than  does  MD, 
some  of  the  advantages  of  which  (e.g.,  the  ability  to  calculate 
transport  properties)  disappear  when  standard  thermostats 
and  barostats  are  invoked  due  to  their  inherently  stochastic 
foundation.51 
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APPENDIX  A:  AVERAGE  MANY-BODY  CONTRIBUTION 
PER  PAIR 

•jj  (2)  _ 

Here  we  prove  the  claim  that  UA/J-  UAH  as  defined  by 
Eqs.  (13)  and  (14)  in  the  text  gives  an  average  many -body 
contribution  per  pair  that,  when  summed  over  all  A  and  B, 
returns  the  total  configurational  energy  of  the  system  under 
the  assumption  that  many-body  effects  terminate  at  third  or¬ 
der.  The  total  energy  as  a  sum  of  two-  and  three-body  inter¬ 
actions  is  given  by 
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where  t/1:'1  represents  the  y-body  contribution  to  the  total 
energy  of  a  v-particle  system.  Let  Ubox.A  equal  the  box  en¬ 
ergy  with  molecule  A  removed. 
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and  assume  an  analogous  relation  for  Uhox_B.  When  both  A 
and  B  are  removed  together, 
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From  Eq.  (A3)  one  can  define  an  effective  pair  energy  in¬ 
cluding  an  average  many-body  correction  as 
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equivalent  to  Eq.  (14)  in  the  text.  Note  that 
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by  comparison  with  Eq.  (Al).  If  the  difference  between  ef¬ 
fective  pair  and  standard  pair  energies  is  defined  as  SUn 
=  U.'.-U]j',  then  the  average  value  of  these  three  quantities 
in  a  A-body  system  is  given  by 
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The  box  energy  expressed  in  terms  of  Eqs.  (A  10)  and  (All) 
is  then 
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By  equating  the  box  energy  of  Eq.  (A5)  with  that  of  Eq. 
(A12)  one  finds  that 
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the  relation  to  be  proved. 


APPENDIX  B:  DERIVATION  OF  EQUATIONS  (43)  and 
(44) 

Treating  (U)  as  a  function  of  /3=(kT)~1  and  applying  the 
chain  rule. 
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dT  T 


(Bl) 


(B2) 


Recall  that  an  ensemble-averaged  quantity  X  is  defined  as 
f  fXewdrdV 


<*>= 


f  JewdrdV 


(B3) 


where  W  represents  the  thermodynamic  weight  appropriate 
to  the  ensemble  being  sampled;  for  the  isothermal-isobaric 
ensemble,  W  is  given  in  Eq.  (16)  of  the  text.  For  the  sake  of 
brevity  we  drop  the  explicit  indication  of  constant  P  in  the 
derivatives  of  Eq.  (Bl)  and  the  Ain  V  factor  appearing  in  the 
exponential  of  Eq.  (16),  then  write  out  explicitly  the  integrals 
appearing  in  the  right-hand  side  of  Eq.  (Bl), 


J 


d(U) 

dp 


(ff(U  +  FV)Ue-^(u+PV)dTdV)(ffe-/3(u+pv)drdV)  (ffUe^(u+PV)dTdV)(ff(U  +  PV)e~0(u+rv)dTdV) 
{jje-P(u+PV)dTdv)  2  +  (Jfe-^u+PV)dTdV) 2  ' 


(B4) 


Using  Eq.  (B3),  Eq.  (B4)  becomes 

~~  =  -  (u(u  +  PV))  +  (U)(U + PV) 
dp 

=  P({U)(V)  -  (uv))  +  (u)2  -  (u2). 


Combining  Eqs.  (B2)  and  (B5)  gives 


d{U) 

dT 


=  £[P«UV)  -  (U)(V))  +  ( U 2)  -  <t/>2]. 


(B6) 


(B5) 


identical  to  Eq.  (43)  in  the  text.  The  derivation  of  Eq.  (44) 
proceeds  in  a  fashion  perfectly  analogous  to  that  detailed 
here  for  Eq.  (40). 
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